
function y_hat = local_linear_rdd(y1,x1,n1,y2,x2,n2,y3,x3,n3,y4,x4,n4,g1,g2,k)

% calculates some numbers once
d = size(x1,2);
ng1 = size(g1,1);
ng2 = size(g2,1);

% standardizes data for each region
std_vcv1=chol(cov(x1));
std_mu1=mean(x1);
std_x1=(x1-repmat(std_mu1,[n1,1]))/std_vcv1;

std_vcv2=chol(cov(x2));
std_mu2=mean(x2);
std_x2=(x2-repmat(std_mu2,[n2,1]))/std_vcv2;

std_vcv3=chol(cov(x3));
std_mu3=mean(x3);
std_x3=(x3-repmat(std_mu3,[n3,1]))/std_vcv3;

std_vcv4=chol(cov(x4));
std_mu4=mean(x4);
std_x4=(x4-repmat(std_mu4,[n4,1]))/std_vcv4;

% calculates bandwidth for each region
bwidth(1)=k*(n1^(-(1/(4+d))));
bwidth(2)=k*(n2^(-(1/(4+d))));
bwidth(3)=k*(n3^(-(1/(4+d))));
bwidth(4)=k*(n4^(-(1/(4+d))));

% bwidth(1)=k*((4/(2*d+1))^(1/(d+4)))*(n1^(-(1/(4+d))));
% bwidth(2)=k*((4/(2*d+1))^(1/(d+4)))*(n2^(-(1/(4+d))));
% bwidth(3)=k*((4/(2*d+1))^(1/(d+4)))*(n3^(-(1/(4+d))));
% bwidth(4)=k*((4/(2*d+1))^(1/(d+4)))*(n4^(-(1/(4+d))));
disp(['bandwidth = ' , num2str(bwidth)]);

% runs LL regression for each grid point
tic;
y_hat=zeros(ng1,ng2);
for i=1:ng1  
    for j=1:ng2
    
        
        if ((g1(i,:)<.5) && (g2(j,:)<.5))
        
            % transforms grid point to standardized scale
            tg=([g1(i,:),g2(j,:)]-std_mu1)/std_vcv1;        
            % generates matrices once
            x_diff=std_x1-repmat(tg,[n1,1]);
            % calculates weights
            t_wgt = mvnpdf((x_diff)/bwidth(1));
            wgt = t_wgt/sum(t_wgt,1);
            % calculates LL coefficients
            z = [ones(n1,1),x_diff];
            zw = z.*repmat(wgt,[1,size(z,2)]);
            b = (zw'*z)*zw'*y1;
            y_hat(i,j) = b(1);
            clear x_diff t_wgt wgt z zw b;
        
        elseif ((g1(i,:)>.5) && (g2(j,:)<.5))

            % transforms grid point to standardized scale
            tg=([g1(i,:),g2(j,:)]-std_mu2)/std_vcv2;        
            % generates matrices once
            x_diff=std_x2-repmat(tg,[n2,1]);
            % calculates weights
            t_wgt = mvnpdf((x_diff)/bwidth(2));
            wgt = t_wgt/sum(t_wgt,1);
            % calculates LL coefficients
            z = [ones(n2,1),x_diff];
            zw = z.*repmat(wgt,[1,size(z,2)]);
            b = (zw'*z)*zw'*y2;
            y_hat(i,j) = b(1);
            clear x_diff t_wgt wgt z zw b;
            
        elseif ((g1(i,:)<.5) && (g2(j,:)>.5))

            % transforms grid point to standardized scale
            tg=([g1(i,:),g2(j,:)]-std_mu3)/std_vcv3;        
            % generates matrices once
            x_diff=std_x3-repmat(tg,[n3,1]);
            % calculates weights
            t_wgt = mvnpdf((x_diff)/bwidth(3));
            wgt = t_wgt/sum(t_wgt,1);
            % calculates LL coefficients
            z = [ones(n3,1),x_diff];
            zw = z.*repmat(wgt,[1,size(z,2)]);
            b = (zw'*z)*zw'*y3;
            y_hat(i,j) = b(1);
            clear x_diff t_wgt wgt z zw b;            
            
        elseif ((g1(i,:)>.5) && (g2(j,:)>.5))

             % transforms grid point to standardized scale
            tg=([g1(i,:),g2(j,:)]-std_mu4)/std_vcv4;        
            % generates matrices once
            x_diff=std_x4-repmat(tg,[n4,1]);
            % calculates weights
            t_wgt = mvnpdf((x_diff)/bwidth(4));
            wgt = t_wgt/sum(t_wgt,1);
            % calculates LL coefficients
            z = [ones(n4,1),x_diff];
            zw = z.*repmat(wgt,[1,size(z,2)]);
            b = (zw'*z)*zw'*y4;
            y_hat(i,j) = b(1);
            clear x_diff t_wgt wgt z zw b;               
            
        end
            
    end 
end
toc;

% saves output
% save output y_hat g1 g2;

% plots the surface
% surf(g1,g2,y_hat');

end






























